Method and apparatus for imaging absorbing objects in a scattering medium

ABSTRACT

A method and processing device are presented for reconstructing an absorption and/or scattering image of a region of interest inside a scattering medium. A mathematical model is provided being representative of a relation between the distribution of the intensity and phase of electromagnetic radiation components scattered from a medium and a certain attenuation factor, which is function of spatial variations of scattering and absorption coefficients of the medium. The mathematical is used for processing a map of distribution of the intensity of electromagnetic radiation components scattered from known locations within the region of interest, thereby producing a halftone pattern of the region of interest.

FIELD OF THE INVENTION

This invention is generally in the field of imaging techniques, andrelates to a method and apparatus for real-time imaging of absorbingobjects within a scattering medium. The invention utilizes ultrasoundtagging of light, and is particularly useful for non-invasivedetection/measurements of absorbing agents, such as hemoglobin, inbiological tissues.

BACKGROUND OF THE INVENTION

In recent years, much effort has been devoted to map the inside ofdiffusive bodies using ultrasound and electromagnetic waves. If a sound(or ultrasound) wave is located inside a scattering medium andconcurrently a continuous electromagnetic wave (such as a laser lightbeam) crosses said medium and is strongly diffused thereby, theelectromagnetic wave frequency is shifted by the sound frequency(acousto-optic effect) at the location of the sound (or ultrasound)wave, while in the other regions, the frequency of the light isunchanged. The electromagnetic waves with the shifted or taggedfrequencies are detected. Since the location of the ultrasound wavesinside the medium, and consequently, the locations of interactionbetween the ultrasound and electromagnetic waves can easily bedetermined, a 3-D representation of the medium can be obtained.

Lev A. et al. “Ultrasound probing of the banana photon distribution inturbid media”, Biomedical Optoacoustics II, San Jose, Calif., USA, 23-24Jan. 2001, vol. 4256, pp. 233-240 discloses the possibility ofultrasound tagging of light to map the photon density inside solidturbine media. The modulation of the optical field transmitted through ascattering medium by an ultrasound beam is also disclosed in M. Kampe etal. “Acousto-optic tomography with multiply scattered light”, OpticalSociety of America, 1997, pp. 1151-1158.

It has been shown [Optics Letters, Lev et al, March 2000] that thetechnique of ultrasound tagging of light provides for locating anelectromagnetic wave absorbing object within a non-absorbing, diffusivemedium. However, it appears that when the are several absorbing objectsor the single object has a pattern of absorbing locations within thediffusive medium, a correlation between the absorption of the differentobjects/locations occurs, and the so-obtained 3-D representation isinsufficient to provide an exact picture of the absorbing pattern withinthe medium, and data reconstruction is thus required.

Image reconstruction techniques typically used with optical measurementsutilize inverse scattering algorithms [S. R. Arridge and J. C. Hebden,“Optical imaging in medicine II. Modeling and reconstruction”, Physicsin Medicine and Biology 42, 841-853 (1997)]. In these methods lightscattered from the medium is detected enabling a two dimensional datarepresentation. This two-dimensional data is then reconstructed into athree-dimensional pattern of absorption (or scattering). The results ofsuch techniques are limited by several factors: the optical measurementmethods are very sensitive to boundary conditions (sensors or sourcespositions), the data reconstruction requires long computation time, andthe image resolution is relatively low, e.g., general not exceeding 5 mmin the case of optical tomography.

SUMMARY OF THE INVENTION

There is accordingly a need in the art to enable 2D or 3D imaging of amabsorbing pattern within a light scattering medium, by providing a novelmethod and apparatus for acquiring data representative of the intensitydistribution of the light response at known locations (voxels) insidethe scattering medium, and processing these data to reconstruct theimage including an absorbing and/or scattering patter in this medium.

The present invention takes advantage of the acoustic tagging of thelight technique enabling identification of the locations inside themedium corresponding to the detected scattered light components. Inother words, irradiation of the medium with acoustic waves is usedsolely for the purpose of mapping the detected scattered lightcomponents. The acoustic tagging of light provides for obtaining the 3-Drepresentation (map) of the medium in real-time. Although thisrepresentation by itself is indicative of the location of an absorbingregion, it does not provide for identifying the absorbing pattern withinthis region (inter alia, because of the cross talk between the differentpaths that light tends to take within the patterned region).

The main idea of the present invention consists of providing an image(halftone pattern) of absorbing and/or scattering inhomogenenitiesinside the scattering medium by processing the map (pattern) of thedetected light signal (e.g., associated with the locations ofinteraction between acoustic and electromagnetic waves). To this end,the map of the detected light signals is considered as a function of theattenuation factor, which, in turn, depends on both the scattering andabsorption properties. In other words, the attenuation factor isrepresentative of the changes in the detected light intensity patternscattered from a scattering medium caused by the absorbing patternthereinside and/or by spatial variations of the scattering properties ofthe medium.

To achieve the above, the measured intensity pattern γ(r₁) for eachposition t inside the medium (e.g., tagged location) is, on the onehand, related to the absorption and scattering coefficients at thatlocation, and, on the other hand, related to the diffusion properties ofthe medium defined by the photons propagation inside the medium. Thesetwo relations provide the basis of the reconstruction algorithmaccording to the invention.

The attenuation factor μ(r₁) is given by:μ=√{square root over (3μ_(a)(μ′_(s)+μ_(a)))}wherein μ_(a) is the absorption coefficient, and μ′_(s) is the reducedscattering coefficient [“Tissue Optics” V. Tuchin, SPIE Press (2000)].

If the scattering coefficient does not change drastically from one pointto another, the determination of the attenuation factor is directlyrelated to the absorption coefficient. In the cases where the scatteringcoefficient varies drastically in space, both coefficients should bedetermined independently.

The present invention also takes into consideration the relativelocations of the light source and detector with respect to each other,as well as the number of sources and detectors used in the measuringdevice.

Thus, according to one aspect of the present invention, there isprovided a method of reconstructing an image of a region of interestinside a scattering medium, the method utilizing a map of distributionof the intensity of electromagnetic radiation components scattered fromtagged locations in said region of interest, said map corresponding tocertain relative positions of at least one source of the electromagneticradiation (14) and at least one detector (16) of the scatteredelectromagnetic radiation with respect to one another and with respectto the tagged locations inside the medium, the method beingcharacterized in that:

-   -   data indicative of said relative positions of the at least one        electromagnetic radiation source and at least one scattered        radiation detector is utilized to provide a mathematical model        that presents a relation between a map of distribution of the        intensity of electromagnetic radiation components γ scattered        from tagged locations in a medium and a halftone map of an        attenuation for μ(r), which is a function of spatial variations        of scattering and absorption coefficients of the medium, said        relation being indicative of a degree of inhomogeneity of at        least one of the scattering and absorption coefficient within        the medium; and,    -   said mathematical model is utilized for processing the map of        the distribution of the electromagnetic radiation components to        thereby obtain a halftone pattern of said region of interest.

According to another aspect of the present invention, there is provideda processing device for use with a measuring apparatus which is operableto detect the intensity of light scattered from tagged locations in ascattering medium, the processing device comprising inter alia a memoryutility for storing a data indicative of relative positions of at leastone light source and at least one detector with respect to one anotherand with respect to the tagged locations inside the medium as used insaid measuring apparatus, and comprising a data processing and analyzingutility for processing input data and generating data indicative of theprocessing results, the processing device being characterized in that:

-   -   said processing device is operable to utilize said data        indicative of the relative positions to create a predetermined        mathematical model and store it in the memory utility, said        mathematical model presenting a relation between a map of the        measured intensities distribution γ and a halftone map of a        certain attenuation factor μ(r), which is a function of the        spatial variation of scattering and absorption coefficients of        the medium, said relation being indicative of a degree of        inhomogeneity of at least one of the scattering and absorption        coefficients within the medium; and

said data processing and analyzing utility, is responsive to input data,coming from the measuring apparatus and being indicative of a map ofdistribution of the intensity of light scattered from the taggedlocations within said region of interest, to apply said mathematicalmodel to said measured data, and thereby obtain a halftone pattern ofthe region of interest.

There is also provided according to yet another aspect of the invention,an apparatus for imaging a region of interest in a scattering medium,the apparatus comprising:

-   -   (i) a measuring apparatus operable to measure electromagnetic        radiation components scattered from different identifiable        locations inside the region of interest, and generate measured        data indicative of a map of the intensity distribution of the        electromagnetic radiation within the region of interest; and    -   (ii) a processing device having a memory unit for storing a        predetermined mathematical model presenting theoretical data        representative of a relation between distribution of the        intensity of electromagnetic radiation scattered from a        scattering medium and a certain attenuation factor, which is a        function of the spatial variations of scattering and absorption        coefficients of the medium; and a processing utility responsive        to input data including said measured data and preprogrammed to        process the measured data with the mathematical model to thereby        obtain a halftone pattern of the region of interest and enable        imaging of the region of interest.

BRIEF DESCRIPTION OF THE DRAWINGS

In order to understand the invention and to see how it may be carriedout in practice, a preferred embodiment will now be described, by way ofnon-limiting example only, with reference to the accompanying drawings,in which:

FIG. 1 is a schematic illustration of an exemplary apparatus constructedand operated for imaging a region of interest by ultrasound tagging oflight, which is suitable to be used for the purposes of the presentinvention;

FIG. 2 more specifically illustrates a processing device according tothe invention used with the apparatus of FIG. 1;

FIGS. 3A and 3B schematically illustrate the conventional imagingtechnique utilizing the ultrasound tagging of light;

FIG. 4 schematically illustrates the cross-talk problem to be solved bythe present invention;

FIG. 5A to 5C schematically illustrates three possible examples of anelectromagnetic source/detector arrangement suitable to be used in thepresent invention;

FIG. 6 schematically illustrates another example of a measurement devicesuitable to be used with the present invention;

FIG. 7 illustrates a flow diagram of the main operational steps of themethod according to the invention; and

FIGS. 8A to 8C illustrate simulation results of applying the method ofthe present invention to obtain a halftone pattern of a region ofinterest in a scattering medium.

DETAILED DESCRIPTION OF THE INVENTION

Referring to FIG. 1, there is illustrated a measurement apparatus 10 forobtaining a map of intensity of the distribution of light scattered froma plurality of known locations with a region of interest (ROI) in ascattering medium (for example, biological tissue) utilizing ultrasoundtagging of light. The apparatus 10 comprises such main constructionalparts as an ultrasound (generally, acoustic) unit 12 coupled to themedium, an illuminator 14 (constituting an electromagnetic radiationsource) optically coupled to the medium, a detector 16, whose output isconnectable to a control unit 18. The apparatus 10 also comprises aphase and frequency control utility 20, which, in the present example,is associated with the ultrasound unit 12, which comprises an acousticor ultrasound generator 22 (possibly including an electronic beamforming unit and an array of amplifiers), and a transducer arrangement24. The operation of the ultrasound unit 12 is aimed at delivering theproper ultrasound wave within the region of interest in the medium.

A function generator 15 operates the ultrasound unit 12 and ananalog-to-digital-converter (card) 19 via a triggering signal TS.Concurrently, the generator 22 transmits an electrical signal to thetransducer arrangement 24 through the phase control utility 18 tothereby actuate one or more transducers to transmit, respectively, oneor more ultrasound signals 25 into a region of interest in the medium.

The illuminator 14 comprises one or several laser devices 26 generatingincident radiation of at least one wavelength (for example, in a rangeof 690-900 nm in order to detect hemoglobin in biological tissue), whichpropagates towards the region of interest. Laser light is diffused(scattered) by the medium, and the diffused light 27 interacts with theultrasound wave 25, and the detector 16 detects the signal resultingfrom this interaction. The electric output of the detector 16 isdirected to the analog-to-digital converter through a band-pass filterand amplifier 28, to thereby produce corresponding digital signalsreceived by the control unit 18.

The control unit 18 comprises a processing devices 30 and 32, and adisplay 34 for displaying an image of the region of interest. Theprocessing device 30 is pre-programmed to process and analyze the datareceived from the detector 16 in the conventional manner, namely, byapplying a power spectrum operation to the measured data and identifyingvariations in light intensity at different frequencies, to therebygenerate measured data indicative thereof. The operational principles ofthe processor 30 do not form part of the present invention, andtherefore need not be described in detail. The processing device 32 isconnectable to the output of the device 30 for further processing themeasured data according to the invention.

As shown in FIG. 2, the processing device 32 is typically a computerdevice comprising a data input utility 36 for receiving the measureddata from the processor 30, a memory utility 38, a data processing andanalyzing utility 40, and a data output utility 42 connectable to thedisplay 34 and/or any other data presenting or data transmittingutility. The memory utility 38 stores a predetermined mathematicalmodel, and the utility 40 is preprogrammed to process the measured datawith the mathematical model, as will be described further below.

It should be understood that input and output data may be of the kind tobe transmitted through wires or wireless to other devices. The input andoutput utilities are constructed accordingly. It should also beunderstood that the processing device 32 can be a separate unitconnectable to the measurement apparatus 10 that includes the controlunit 18 with the processor 30 as its constructional part.

Turning back to FIG. 1, it should be noted that the phase controlutility 20 may be alternatively associated with the illuminator 14, andmay be a part of the function generator 15. In this case, the functiongenerator operates to modulate the output intensity of the laser using aphase modulation scheme. The measurement device may be of any knownkind, for example as disclosed in U.S. Pat. No. 6,041,248, and isdisclosed in the co-pending application assigned to the assignee of thepresent application. The construction and operation of the measurementdevice do not form part of the present invention, and therefore are notspecifically described, except to note the following:

The interaction between the light wave and the ultrasound results inthat the frequency/phase of light is shifted by the frequency/phase ofthe ultrasound, and the presence of an absorbing agent in the scatteringmedium can be determined from the change in the intensity distributionof the frequency/phase shifted light signals. The light source, theprobed region, and the detector do not have to be specifically alignedwith each other, and can have any geometric configuration, provided thatenough photons reach the detector. This allows multiple-source/detectorconfigurations, with an increase in the signal to noise ratio and betterlight filling of the tissues.

The interaction is as follows: The light source emits light of frequencyω into the probed region (region of interest). The ultrasound pulses offrequency ΩUS are transmitted into the probed region. The currentlocation(s) of the interaction in the X-Y plane is defined by thecurrent location of the transducer(s), and along the Z-axis by the phaseof ultrasound pulses. The ultrasound modulated light having a shiftedfrequency ω+ΩUS, and non-modulated light having the frequency ω arereceived by the detector, which mixes them and generates a signalmodulated at the ultrasound frequency. Data indicative of the detectedsignals is processed in the device 30 to obtain the measured data in theform of a map of the modulated light intensity distribution within theregion of interest.

It should be noted that the present invention is not limited to theabove-described technique of obtaining data representation (e.g., 3Drepresentation) of the scattered light map. Any other suitable techniqueand device can be used for obtaining such measured data to be furtherprocessed by the present invention. For example, the technique of theabove-indicated U.S. Pat. No. 6,041,248 is suitable.

FIG. 3A schematically illustrates a region of interest (ROI), which is ascattering medium (e.g., biological tissue), and which has an absorbingregion (AR) thereinside. FIG. 3B schematically illustrates the measureddata (MD) obtained with the processor 30. The measured data is the lightintensity map typically shaped like a banana between the source and thedetector [S. C. Feng, F. Zeng and B. Chance, SPIE Vol. 2389,pp 54-63.1995], and is distorted by the presence of the absorbing region. It isthus evident that in order to obtain a halftone pattern (image) of theabsorbing region, further processing of the measured data is needed.

The ultrasound tagging of light allows to obtain a data representationof the region of interest (e.g., three-dimensional representation).However, these data do not give rise directly to the absorptiondistribution (pattern) because of the cross talk between the differentpaths that photons can take. This is illustrated in FIG. 4, showing aregion of interest 42 in a scattering medium having two spaced-apartabsorbers 44 a and 44 b. Photons generated by a light source 46propagate through the region of interest in two photon paths towards adetector 48. Although a separate tagging zone is located at the absorber44 b, the photons crossing the absorber 44 a will influence the taggedsignal of the absorber 44 b at the detector. Thus there is a need for areconstruction method that will enable to eliminate these cross-talkeffects.

According to the invention, the data reconstruction method is based onconsidering the measured intensity pattern γ(r_(t)) for each position tof the tagging location, as a function of the attenuation factorμ(r_(t)). The latter depends on both the scattering and absorptionproperties of the region of interest, and is actually representative ofthe changes in the detected light intensity pattern caused by theabsorbing pattern in the region of interest. Hence, the measuredintensity pattern γ(r_(t)) for each tagged location t is, on the onehand, related to the absorption coefficient at that location, and, onthe other hand, related to the diffusion properties of the mediumdefined by the light propagation inside the medium.

The attenuation factor is determined as follows:μ=√{square root over (3μ_(a)(μ_(s)+μ_(a)))}wherein μ_(a) is the absorption coefficient and μ′_(s) is the reducedscattering coefficient.

If the scattering coefficient does not change drastically from onelocation to another, the determination of the attenuation factor isdirectly related to the determination of the absorption coefficient. Inthe medium of the kind where the scattering coefficient variesdrastically in space, both the absorption and scattering coefficientsshould be determined independently. In the case of brain trauma orstroke, or in the case of cancer detection based on the angiogenesisprocess, the effect of variations in the absorption pattern of themedium due to the presence of blood vessels is much more dominant thanthat of the variations in the scattering properties. Therefore, in thatcase, the effect of changes in scattering properties within the mediumcan generally be neglected. As for effects in bones for examples, thescattering coefficient's variations between the locations of thecortical bone and the trabecular bone, or between the locations of thesane and diseased bones (e.g., in the case of osetoporosis), are quiteimportant. Accordingly, in these cases, changes in scattering propertiesshould also be taken into consideration.

It is known [A. Ishimaru, “Wave Propagation and Scattering in RandomMedia”, Academic Press, (1978)] that electromagnetic waves' propagationin scattering media can be considered (in a first but well-establishedapproximation) as experiencing a diffusion process (the so-called P₁approximation). Therefore, a diffusion equation for the photons' number(or equivalently, the photon or energy density) considering theattenuation factor along three orthogonal axes is as follows:

 (∇²−μ²(r))U(r)=S(r)  (1)

wherein U(r) is the intensity of the scattered light (electromagneticwave); S(r) characterizes the sources distribution (i.e., thedistribution of the illuminating intensity); and μ(r) is the attenuationfactor or the so-called extinction coefficient.

As indicated above, when the light with a frequency ω interacts with theultrasound (acoustic) beam of a frequency Ω, a new frequency ω+Ω iscreated (with a certain tagging efficiency η). This, however, does notaffect the scattering properties of the medium, and thus the diffusionequation (1) remains the same.

Turning back to FIG. 4, let us consider the following geometricconfiguration: the light source 46 is located at the surface locationr_(s), the detector 48 is located at r_(d), and the ultrasound beam islocated at r_(t) inside the medium. Here, the index “t” denotes thetagging zone, namely, the location in the region of interest where lightfrequency is modulated.

The probability for a photon to migrate from the first location r₁ tothe second location r₂ is denoted by Γ(r₁,r₂). The probability to detecta modulated photon is proportional to the product of the probability tomigrate from the source to the tagging zone, the tagging efficiency(defined as the percentage of photons crossing the ultrasound wave thateffectively get tagged), and the probability to migrate from the taggingzone to the detector. The tagging efficiency is defined as thepercentage of photons crossing the ultrasound wave that effectively gettagged

Therefore, the detected light intensity distribution γ(r_(t)) is asfollows:γ(r_(t))≅Γ(r_(s),r_(t))Γ(r_(t),r_(d))ηS₀  (2)wherein S₀ is the source intensity (i.e., the intensity of the incidentlight).

In the above equation, intensity of photons that several timesinteracted with the acoustic waves has been neglected, since theproportion scales of these photons is small.

In general, when the detector, the source and the tagging zone are alllocalized in three different points, the probabilities Γ are reduced tothe Green function G(r) of the equation (1), that isγ(r_(t))≅(2π)²G(r_(s)−r_(t))G(r_(t)−r_(d))ηS₀  (3)

In the general case, where the light source is not a point source and/orthe detector is not a point detector, the above equation (3) must bereplaced by the following:γ(r _(t))=qj ₀ ∫∫dA _(s) dA _(d)Γ(r _(s) ,r _(t))Γ(r _(t) ,r _(d))=ηS₀Γ_(s)(r _(t))Γ_(d)(r _(t))  (4)wherein A_(s) and A_(d) represent, respectively, the source and detectorarea, and Γ_(s)(r_(t)) and Γ_(d)(r_(t)) are the probabilities of thephoton to migrate from, respectively, the source to the tagging locationr_(t) and from the tagging location r_(t) to the detector. It should benoted that since Γ(r₁, r₂)=Γ(r₂,r₁), the probability Γ_(d)(r_(t)) couldalso be referred to as the probability of the photon to migrate from thedetector to the tagging location r_(t).

Accordingly, we have: $\begin{matrix}{\frac{\nabla^{2}\sqrt{\gamma}}{\sqrt{\gamma}} \equiv {{- {\frac{1}{4}\left\lbrack {\frac{\nabla\Gamma_{s}}{\Gamma_{s}} - \frac{\nabla\Gamma_{d}}{\Gamma_{d}}} \right\rbrack}^{2}} + {\frac{1}{2}\left\lbrack {\frac{\nabla^{2}\Gamma_{s}}{\Gamma_{s}} + \frac{\nabla^{2}\Gamma_{d}}{\Gamma_{d}}} \right\rbrack}}} & (5)\end{matrix}$

The above equation (5) is an exact algebraic relation without any hiddenphysical assumption, and is therefore always valid.

As there is no light source inside the medium, both probabilitiesΓ_(s)(r_(t)) and Γ_(d)(r_(t)) are solutions of the diffusion equation(1) for any location r_(t) inside the medium, but without thesource-associated term:∇²Γ(r _(t))=μ²(r _(t))Γ(r _(t))  (6)

Hence, equation (5) can be rewritten as follows: $\begin{matrix}{\frac{\nabla^{2}\sqrt{\gamma}}{\sqrt{\gamma}} = {{- {\frac{1}{4}\left\lbrack {\frac{\nabla\Gamma_{s}}{\Gamma_{s}} - \frac{\nabla\Gamma_{d}}{\Gamma_{d}}} \right\rbrack}^{2}} + {\mu^{2}\left( r_{t} \right)}}} & (7)\end{matrix}$wherein all the derivatives are taken with respect to the taggedlocations r_(t). From this point of view, the only approximation thatleads to equation (7) is the equation (2), i.e., the multiplicationassumption.

Reference is now made to FIGS. 5A-5C, illustrating three possibleexamples of an electromagnetic source/detector arrangement.

In the example of FIG. 5A, the source 46A and the detector 48A are atthe same location and have the same dimensions. In this case, equation(7) can be simplified and yields the following: $\begin{matrix}{\frac{\nabla^{2}\sqrt{\gamma\left( r_{t} \right)}}{\sqrt{\gamma\left( r_{t} \right)}} = {\mu^{2}\left( r_{t} \right)}} & (8)\end{matrix}$

In the example of FIG. 5B, the source and the detector 46B and 48B arespaced from each other, and a certain angle θ therefore exists betweenthe vectors (r_(t)−r_(s)) and (r_(t)−r_(d)) presenting the central axesof incident and scattered light components propagating to and from theabsorbing tagged region AR. In this case, the second term in theequation (5) above cannot be neglected any more and must be estimated.In the second order approximation, we can evaluate Γ as if the medium ishomogenous. In that case, the analytical expression of Γ is well-known[A. Ishimaru, “Wave Propagation and Scattering in Random Media”,Academic Press, (1978)] and therefore the expression$\frac{\nabla\Gamma}{\Gamma}$can be developed in powers of 1/r. For distances r which are larger than1/{overscore (μ)}, {overscore (μ)} being the average value of theattenuation factor μ(r), we have the general expression: $\begin{matrix}{{{\frac{\nabla\Gamma}{\Gamma}(r)} \sim {{{- {\mu(r)}}\frac{r}{r}} - \frac{r}{{r}^{2}}}}{{and}\quad{therefore}\text{:}}\frac{\nabla^{2}\sqrt{\gamma}}{\sqrt{\gamma}} \cong {{\mu^{2}\left( r_{t} \right)}\quad\frac{\left\lbrack {1 + {\cos\quad\theta}} \right\rbrack}{2}}} & (9)\end{matrix}$It should be understood that in equation (9), r represents any positionin the medium, and therefore states a physical approximation that is notlinked with the ultrasound tagging of light. In the last equation, thevalue of r is replaced by r_(t) thereby taking into account theultrasound tagging.

In the case of spaced-apart source and detector, the image obtained atthe detector should not be drastically different from that obtained withthe source and detector located adjacent to each other. Therefore,equation (8) should just be corrected as presented in the last equation.This correction is estimated by taking the solution of the problem as ifthe medium is homogenous, and is added to equation (8).

Therefore, equation (8) can be generalized to the case where the sourceand the detector are at different places, as follows: $\begin{matrix}{{\mu^{2}\left( r_{t} \right)} = {\frac{2}{1 + {\cos\quad\theta}}\frac{\nabla^{2}\sqrt{\gamma}}{\sqrt{\gamma}}}} & (10)\end{matrix}$

Equation (10) returns to equation (8) when θ→0, namely when the distancebetween the source and the detector is much smaller than the distancebetween the source/detector and the tagging zone (the acoustic wave).

FIG. 5C exemplifies the situation with an array of sources 46C₁-46C₅(generally at least two) and an array of detectors 48C₁-48C₅ (similarly,at least two detectors). This case is the most general, and requires aniterative algorithm consisting of estimating the light intensitydetected by each detector considering that this detected light intensityis produced by scattered light components of all the sources. Thus, theiterative algorithm includes the following steps:

Step 1: Equivalent location of the source r_(s) is evaluated as follows:$\begin{matrix}{r_{s} = \frac{\sum\limits_{j = 1}^{N_{s}}\quad{S_{j}r_{j}}}{\sum\limits_{j = 1}^{N_{s}}\quad S_{j}}} & (11)\end{matrix}$wherein S_(j) is the intensity of the j^(th) source at the locationr_(j).

From the above-described example of FIG. 5B, μ^((i))(r_(t)) being thefirst approximation of μ(r_(t)) can be evaluated using the followingequation: $\begin{matrix}{{{+ 2}d\quad{\mu^{(1)}\left( r_{t} \right)}} = {\frac{1}{N_{d}}{\sum\limits_{i = 1}^{N_{d}}\quad{\frac{2}{1 + {\cos\quad\theta_{i}}}\frac{\nabla^{2}\sqrt{\gamma_{(i)}\left( r_{t} \right)}}{\sqrt{\gamma_{(i)}\left( r_{t} \right)}}}}}} & (12)\end{matrix}$wherein N_(d) is the total number of detectors, γ_((i))(r_(t)) is thelight intensity returned from the tagging location r_(t) to the i^(th),and θ_(i) is the angle between the vectors r_(t)−r_(i) and r_(t)−r_(s).

Step 2: The diffusion equation (1) is solved using μ⁽¹⁾(r_(t)) insteadof μ(r_(t)). To this end, two equations must be solved: one for thesources and one for the detectors. The boundary conditions must bechosen so as to match the sources and the detectors' arrangement:−∇²Γ_(s) ⁽¹⁾(r)+μ₍₁₎ ²(r)Γ_(s) ⁽¹⁾(r)=S(r)  (13)−∇²Γ_(d) ⁽¹⁾(r)+μ₍₁₎ ²(r)Γ_(d) ⁽¹⁾(r)=S(r)  (14)

Step 3: Now, the new value for μ(r_(t)), i.e., μ₍₂₎(r_(t))is obtainedusing the following: $\begin{matrix}{{\mu_{(2)}^{2}\left( r_{t} \right)} = {{\mu_{(1)}^{2}\left( r_{t} \right)} + {\frac{1}{4}\left\lbrack {\frac{\nabla\Gamma_{s}^{(1)}}{\Gamma_{s}^{(1)}} - \frac{\nabla\Gamma_{d}^{(1)}}{\Gamma_{d}^{(1)}}} \right\rbrack}^{2}}} & (15)\end{matrix}$

Step 4: Steps 2 and 3 are repeated until a precise enough value forμ(r_(t)) is obtained. The repetitive process can for example proceeduntil the difference in the value of μ(r_(t)) from two consecutiveiterations is smaller than ε (in absolute value).

The method of the present invention can be applied directly to thedetermination of biological tissues' saturation. By choosing twodifferent wavelengths of incident light in the near-infrared spectralregion (670 nm to 900 nm), it is possible to obtain the value ofμ(r,λ)for both wavelengths. The absorption of such a medium isessentially due to oxyhemoglobin and deoxyhemoglobin, and the manner ofextracting the oxygen saturation of tissues from the measured scatteredlight is known in the art. If scattering properties of the medium arenot varying strongly from one place to another, the ratio ofμ(r,λ₁)/μ(r,λ₂) is directly proportional to the ratioμ_(a)(r,λ₁)/μ_(a)(r,λ₂) [S. J. Matcher, P. Kirpatrick, K. Nahid, M. Copeand D. T. Delpy, SPIE Vol. 2389. pp 486-495, 1995]. This ratio leadsdirectly to the saturation.

In some cases, the knowledge of the attenuation factor is insufficientfor imaging purposes. As indicated above, such a situation occurs with amedium of the kind where the scattering coefficient varies drasticallyin space. In this case, both the absorption and the reduced scatteringcoefficients μ_(a) and μ′_(s) should be determined independently. Morespecifically, the spatial distribution of the scattering coefficientwithin the region of interest is to be determined.

The above can be implemented by using a measurement device somewhatdifferent from the device of FIG. 1. This measurement device isschematically illustrated in FIG. 6, being generally denoted 100. Tofacilitate understanding, the same reference numbers are used foridentifying those components which are common in the devices 10 and 100.Thus, in the device 100, there is additionally provided a functiongenerator 35 that modulates the detector response and provides asynchronized reference, an amplifier 36 to amplify the signal to thedetector, a phase detection system (e.g., a phase-lock loop or PLL) 40and a frequency mixer system 42, which are connected to the detector andto the filter 28. It should be understood that the operations of thefunction generators 15 and 35 are appropriately synchronized. Theoperation of the measurement device 100 is based on irradiating theregion of interest with a high-frequency modulated light, and detectingthe phase variations in the ultrasound tagged light components, i.e.,light components scattered from different locations of the region ofinterest. Thus, in this configuration, the detected light signals comingfrom different spatial coordinates, in addition to being frequencytagged due to the interaction with ultrasound, have certain phase tags,thereby enabling separate identification of the light intensity changesdue to the changes in the scattering properties only.

More specifically, the laser is modulated at a high frequency (typicallybetween 50 MHz and 1 GHz) at the modulation ω. This creates a so-calleddiffusive wave within the scattering medium [“Tissue Optics”, V. Tuchin,SPIE Press (2000)]. The photomultiplier dynodes (detector) are alsomodulated at the frequency ω+δωl+Ω, wherein Ω is the ultrasoundfrequency. At the frequency δωl, the detected signal depends on theexact position of the ultrasound pulse inside the medium. The ultrasoundpulse position determines the distance that the diffusive wavepropagates, and therefore its phase. This phase is determined bycomparing the signal at the frequency δωl with a reference signal at thesame frequency (derived from the same clock as the generator thatmodulates the dynodes). It should be noted that the frequency ω iseither fixed or time dependent. This is important in order to enableobtaining the relation between the ultrasound pulse position and thesignal. It should also be understood that other frequency combinationscould be used in order to obtain a similar result.

The knowledge of both the amplitude and phase of the detected signalgives the knowledge of the complex intensity of the signal. Theamplitude can be determined by any known technique (e.g., U.S. Pat. No.6,041,248) or by the technique of the above-indicated co-pendingapplication utilizing measurement of the power spectrum. The timedependent equation is as follows: $\begin{matrix}{{{- {\nabla^{2}{U(r)}}} + {{\mu^{2}(r)}{U(r)}} - {\frac{1}{D}\frac{\partial U}{\partial t}}} = \frac{S(r)}{D}} & (16)\end{matrix}$wherein $D = \frac{c}{3\left( {\mu_{a} + \mu_{s}^{\prime}} \right)}$is the diffusion coefficient and c is the speed of light in vacuum.

The light source S is now considered as a modulated light source, ratherthan a static light source, as in the previously described examples, andtherefore we have:S(r,t)=S ₀(r)exp(−i(ωt+φ))  (17)

The general form of the electromagnetic wave intensity is as follows:

 U(r,t)=U ₀(r,t)exp(−iωt)  (18)

In that case the diffusion equation becomes: $\begin{matrix}{{\left( {\nabla^{2}{- k^{2}}} \right){U_{0}\left( {r,t} \right)}} = {{- \frac{S_{0}(r)}{D}}\quad{\exp\left( {{- {\mathbb{i}}}\quad\varphi} \right)}}} & (19)\end{matrix}$wherein the complex attenuation coefficient k is as follows:$\begin{matrix}{k^{2} = {\mu^{2} + {{\mathbb{i}}\quad\frac{\omega}{D}}}} & (20)\end{matrix}$

The above equation is conceptually similar to equation (2), anddistinguished from equation (2) in the following: coefficient k replacescoefficient μ; all the values are complex, and therefore, not only theamplitude but also the phase of the detected signal φ is to be measured.

Once k is determined experimentally (both k_(r) and k_(i)), equation(20) is a simple system of two equations with two unknowns: μ_(a) andμ_(s). The solution is given by: $\begin{matrix}\begin{matrix}{\mu_{a} = {\frac{3\quad\omega}{c}\frac{k_{r}^{2} - k_{i}^{2}}{2\quad k_{r}k_{i}}}} \\{\mu_{s}^{\prime} = {{\frac{c}{3\quad\omega}\quad 2\quad k_{r}k_{i}} - {\frac{3\quad\omega}{c}\frac{k_{r}^{2} - k_{i}^{2}}{2\quad k_{r}k_{i}}}}}\end{matrix} & (21)\end{matrix}$

Having scanned the region of interest with the electromagnetic andacoustic waved and measured the intensity distribution of the scatteredlight, either a specific one of the above algorithms or the most generalone is applied to the measured data to obtain a halftone pattern of theregion of interest. If the most general algorithm is used, the real andimaginary parts of the coefficient k² can be retrieved, leading to thetwo equations, which determine the values of μ′_(s) and μ_(a).

Thus, as shown in FIG. 7, to enable the imaging technique of the presentinvention, the measured data indicative of the distribution of theintensity of light scattered from the region of interest is provided. Tothis end, if the spatial distribution of the scattering coefficient inthe region of interest is to be separately determined, theabove-described technique of detecting scattered light using phasemodulated electromagnetic radiation is employed. Then the measured datais processed by applying thereto a selected model (algorithm) inaccordance with the source/detector arrangement (their number andrespective positions), thereby determining the halftone pattern of theregion of interest.

FIGS. 8A-8C illustrate the light tomography simulation results showingthe advantageous features of the present invention. For simplicity, butwithout loss of generality, a 2D medium imaging is simulated. The mediumshown in FIG. 8A is a two-dimensional slightly absorbing but scatteringmedium, and includes an absorbing pattern—the word “SOREQ” which is tobe imaged. Simulated measurement data shown in FIG. 8B, is thenprocessed by the technique of the present invention resulting in thehalftone image of the absorbing pattern as shown in FIG. 8C. In thissimulation, the above described case of the source and detector at thesame location was used.

The simulation is performed as follow: light is injected in the medium(coordinate (16,0) in the picture) and the resulting image is simulated.In the simulation, the photons number in every volume element around rin the 2D medium is split as follows: 1−exp[−μ_(a)], is absorbed by themedium, and the rest is scattered evenly among its four nearestneighbors (rectangular grid). It should be noted that here μ_(a) is notexactly the absorption coefficient but is proportional to the absorptioncoefficient. The tagging process is simulated in the following manner:each photon that crosses the tagging zone is differentiated in thesimulation and is counted separately when it reaches the detector.

It should be noted that the technique of the present invention does notrequire a three-dimensional acoustic scan of the entire body partcontaining the region of interest. It is possible to scan only aselected portion of the body, for example with high resolution (“zoommode”).

It should also be noted that by properly choosing the positions of thesources and detectors, it is possible to minimize the residual term${\frac{1}{4}\left\lbrack {\frac{\nabla\Gamma_{s}^{(1)}}{\Gamma_{s}^{(1)}} - \frac{\nabla\Gamma_{d}^{(1)}}{\Gamma_{d}^{(1)}}} \right\rbrack}^{2}.$This occurs for example when the sources and detectors are symmetricallyarranged, for example, one-detector is located at the center of apolygon composed of sources.

The reconstruction technique of the present invention can be applied incases that are different from imaging of biological tissues. Thefollowing two examples demonstrate the wide range of applications ofthis technique.

In the field of microscopy, for example aimed at examining multilayerprinted circuit boards, high frequency (10 to 200 MHz) ultrasoundtransducers can be used so that the resolution is in the range of tensto hundreds of microns. By choosing the proper light wavelengths,different sorts of dielectric materials (i.e., having differentabsorbing patterns) can be identified in the different layers of thePrinted Circuit Board, although their determination by ultrasonicmethods only might be difficult. Such a microscopy can also be used fordetecting melanoma and other skin lesions by analyzing the precise shapeof the lesion in depth. It is an alternative to confocal microscopy oroptical coherence tomography.

In a different context, the technique can be used for analyzing eitherthe quality or the homogeneity of diffusive bulk-like materials such asfood product or plastics, within their containers. To this end, laserlight is coupled to the container (or directly to the material underevaluation), and an ultrasound pulse is scanned over the material. Inthe case of large containers (centimeters to meters), low frequencyultrasound can be used, which while reducing the resolution, increasesthe penetration depth and the scan speed. Having performed thereconstruction process of the present invention, a three-dimensionalimage of the absorption and/or scattering pattern can be obtained.

Those skilled in the art will readily appreciate that variousmodifications and changes can be applied to the embodiments of theinvention as hereinbefore exemplified without departing from its scopedefined in and by the appended claims.

1. A method of reconstructing an image of a region of interest inside ascattering medium, the method utilizing a map of distribution of theintensity of electromagnetic radiation components scattered from taggedlocations within said region of interest, said map corresponding tocertain relative positions of at least one source of the electromagneticradiation (14) and at least one detector (16) of the scatteredelectromagnetic radiation with respect to one another and with respectto the tagged locations inside the medium the method being characterizedin that: data indicative of said relative positions of the at least oneelectromagnetic radiation source and that at least one scatteredradiation detector is utilized to provide a mathematical model thatpresents a relation between a map of distribution of the intensity ofelectromagnetic radiation components γ scattered from tagged locationsin a medium and a halftone map of an attenuation factor μ(r), which is afunction of spatial variations of scattering and absorption coefficientsof the mediums said relation being indicative of a degree ofinhomogeneity of at least one of the scattering and absorptioncoefficients within the medium; and, said mathematical model is utilizedfor processing the map of the distribution of the electromagneticradiation components to thereby obtain a halftone pattern of said regionof interest.
 2. The method according to claim 1, wherein said relationbetween the intensity distribution and the attenuation factor is basedon the relative position of the at least one source of theelectromagnetic radiation, the at least one detector and at least oneacoustic radiation source used for tagging said locations in the mediumas locations of interactions between the acoustic and electromagneticradiation inside the medium.
 3. The method according to claim 1, whereinsaid processing comprises determining changes in the intensitydistribution of the electromagnetic radiation caused by an absorbingpattern in the region of interest, based on that said intensity for eachlocation in the region of interest is related to the absorptioncoefficient and to the scattering properties of the medium defined bythe electromagnetic propagation inside the medium.
 4. The methodaccording to claim 1, wherein said relation is determined as follows:$\frac{\nabla^{2}\sqrt{\gamma}}{\sqrt{\gamma}} = {{- {\frac{1}{4}\left\lbrack {\frac{\nabla\Gamma_{s}}{\Gamma_{s}} - \frac{\nabla\Gamma_{d}}{\Gamma_{d}}} \right\rbrack}^{2}} + {\mu^{2}\left( r_{t} \right)}}$wherein γ is the detected distribution of the intensity of the,electromagnetic radiation scattered from the tagged locations in themedium; Γ_(s) and Γ_(d) are probabilities of a photon to migrate from,respectively, the electromagnetic radiation source to the taggedlocation inside the medium and from the tagged location in the medium tothe detector; and a μ(r_(t)) is the attenuation factor taken withrespect to tagged locations r_(t) in the medium.
 5. The method accordingto claim 4, wherein the electromagnetic radiation source and thedetector are at the same location, said relation being determined asfollows:$\frac{\nabla^{2}\sqrt{\gamma\left( r_{t} \right)}}{\sqrt{\gamma\left( r_{t} \right)}} = {{\mu^{2}\left( r_{t} \right)}.}$6. The method according to claim 4, wherein the electromagnetic sourceand the detector are spaced from each other, and a certain angle θtherefore exists between central axes of incident and scattered lightcomponents propagating to and from an absorbing region in the medium,said relation being determined as follows:${\mu^{2}\left( r_{t} \right)} = {\frac{2}{1 + {\cos\quad\theta}}{\frac{\nabla^{2}\sqrt{\gamma}}{\sqrt{\gamma}}.}}$7. The method according to claim 4, wherein an array of theelectromagnetic sources and an array of the detectors are used forobtaining the intensity distribution map, said relation being determinedby an iteractive algorithm comprising estimating detection by each ofthe detectors the scattering of the electromagnetic radiation producedby all the sources.
 8. The method according to claim 1, comprisingapplying measurements to the medium to provide said the map ofdistribution of the intensity of electromagnetic radiation componentsscattered from the tagged locations within said region of interest. 9.The method according to claim 8, wherein the measurements are based ontagging the locations within said region of interest irradiated with theelectromagnetic radiation.
 10. The method according to claim 9, whereinsad tagging is based on interaction between the electromagneticradiation with acoustic radiation inside the medium.
 11. The methodaccording to claim 10, wherein said tagging is applied to the amplitudeof the scattered electromagnetic radiation component.
 12. The methodaccording to claim 10, wherein said tagging is applied to the amplitudeand phase of the scattered electromagnetic radiation components.
 13. Themethod according to claim 10, wherein said tagging comprises irradiatingthe region of interest with acoustic and electromagnetic waves toprovide a plurality of locations inside the region of interest where theacoustic and electromagnetic waves interactive with each other, therebyaffecting the parameters of the electromagnetic waves scattered from theregion of interest in accordance with those of the acoustic waves,thereby enabling detection of the scattered electromagnetic waves withthe affected parameters.
 14. A processing device (32) for use with ameasuring apparatus (10), which is operable to detect the intensity oflight scattered from tagged locations in a scattering medium, theprocessing device comprising inter alia a memory utility (38) forstoring data indicative of relative positions of at least one lightsource (14) and at least one detector (16) with respect to one anotherand with respect to the tagged locations inside the medium as used insaid measuring apparatus, and comprising a data processing and analyzingutility (40) for processing input data and generating data indicative ofthe processing results, the processing device (32) being characterizedin that: said processing device is operable to utilize said dataindicative of the relative positions to create a predeterminedmathematical model and store it in the memory utility, said mathematicalmodel presenting a relation between a map of the measured intensitiesdistribution γ and a halftone map of a certain absorption factor μ(r),which is a function of the spatial variation of scattering andabsorption coefficients of the medium, said relation being indicative ofa degree of inhomogeneity of at least one of the scattering andabsorption coefficients within the medium; and said data processing andanalyzing utility is responsive to input data coming from the measuringapparatus and being indicative of a map of distribution of the intensityof light scattered from the tagged locations within said region ofinterest, to apply said mathematical model to said measured data, andthereby obtain a halftone pattern of the region of interest.
 15. Thedevice according to claim 14, wherein said processing and analyzingutility operates to determine changes in the intensity distribution ofthe scattered light caused by an absorbing pattern in the region ofinterest, based on that said intensity for each the locations in theregion of interest is related to the absorption coefficient and to thescattering properties of the medium defined by the light propagationinside the medium.
 16. The device according to claim 14, wherein saidrelation between the intensity distribution and the attenuation factoris related to the relative position of the at least one light source, atleast one detector for detecting the scattered radiation components, andat least one acoustic unit used for tagging said locations in the mediumas locations of interactions between light and acoustic radiation in themedium.
 17. The device according to claim 14, wherein said relation isdetermined as follows:$\frac{\nabla^{2}\sqrt{\gamma}}{\sqrt{\gamma}} = {{- {\frac{1}{4}\left\lbrack {\frac{\nabla\Gamma_{s}}{\Gamma_{s}} - \frac{\nabla\Gamma_{d}}{\Gamma_{d}}} \right\rbrack}^{2}} + {\mu^{2}\left( r_{t} \right)}}$wherein γ is the detected distribution of the intensity of lightscattered from the tagged locations in the medium; Γ_(s) and Γ_(d) areprobabilities of a photon to migrate from, respectively, theelectromagnetic radiation source to the tagged location inside themedium and from the tagged location in the medium to the detector; andμ(r_(t)) is the attenuation factor taken with respect to taggedlocations r_(t) in the medium.
 18. The device according to claim 17,wherein the light source and the detector are at the same location, saidrelation being determined as follows:$\frac{\nabla^{2}\sqrt{\gamma\left( r_{t} \right)}}{\sqrt{\gamma\left( r_{t} \right)}} = {{\mu^{2}\left( r_{t} \right)}.}$19. The device according to claim 17, wherein the light source and thedetector are spaced from each other, and a certain angle θ thereforeexists between central axes of incident and scattered light componentspropagating to and from an absorbing region in the medium, said relationbeing determined as follows:${\mu^{2}\left( r_{t} \right)} = {\frac{2}{1 + {\cos\quad\theta}}{\frac{\nabla^{2}\sqrt{\gamma}}{\sqrt{\gamma}}.}}$20. The device according to claim 17, wherein an array of the lightsources and an array of the detectors are used for obtaining theintensity distribution map, said relation being determined by aniteractive algorithm comprising estimating detection by each of thedetectors the light scattering produced by all the sources.
 21. Thedevice according to claim 14, and also comprising a processorconnectable to the output of said meaning apparatus for receiving andanalyzing data indicative of light components scattered from the mediumto generate said measured data indicative of the map of distribution ofthe intensity of light components scattered from the tagged locations inthe medium.
 22. The device according to claim 14, being connectable to aprocessor associated with said measuring apparatus for receiving andanalyzing data indicative of light components scattered from the regionof interest to generate said measured data indicative of the map ofdistribution of the intensity of light components scattered from thetagged locations in the medium.
 23. The device according to claim 14 andalso comprising a measuring apparatus operable to measureelectromagnetic radiation components scattered from tagged locationsinside the region of interest, and generate measured data indicative ofa map of distribution of the intensity of light scattered from thetagged locations within the region of interest.
 24. The apparatusaccording to claim 23, wherein said measuring apparatus comprises: (i)an acoustic unit operable to transmit acoustic radiation to a pluralityof locations in the medium, and a light source operable to illuminatesaid medium with incident light, to thereby produce light signals, beingmodulated in accordance with the acoustic radiation and allow detectionof light components scattered from tagged locations of interactionsbetween light and acoustic radiation inside the medium; and (ii) adetector unit operable to detect said modulated light signals andgenerate measured data indicative thereof, said measured data being inthe form of a map of intensity distribution of light scattered from thelocations inside the region of interest.
 25. The apparatus according toclaim 24, wherein the light source generates light of one or morewavelengths.
 26. The apparatus according to claim 25, wherein said lightsource produces high-frequency modulated light components and saiddetector is modulated at a nearby high-frequency.
 27. The apparatusaccording to claim 26, wherein said detector is associated with a phasesystem for detecting said modulated light sign and phase variations ofthe light components scattered from different locations in the region ofinterest, said map being indicative of the distribution of thescattering attenuation or within the locations inside the region ofinterest.
 28. The apparatus according to claim 24, and also comprising aphase utility associated with the acoustic unit, thereby allowing forimaging of the region of interest along an axis of transmission of saidacoustic radiation.